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Abstract: 

Hotspots are associated with long- wavelength geoid highs, an association that is 
even stronger when the geoid highs associated with subduction zones are removed. 
We quantify these associations by expanding the hotspot distribution in spherical har- 
monics and calculating correlation coefficients as a function of harmonic degree. The 
hotspot distribution spectrum is essentially white, with peaks at degrees 2 and 6. It 
is correlated positively with the slab residual geoid for degrees 2-6, with low seismic 
velocity in the lower mantle at degree 2, and with low seismic velocity in the upper 
mantle at degree 6. We test a variety of fluid mechanical models for hotspots, includ- 
ing lithospheric delamination and hot plumes, by calculating their predicted dynamic 
geoid responses and comparing them to the observations. These models include the 
effects of temperature dependent rheology. Our preferred hotspot model, based on 
observations of the geoid and seismic tomography, has plumes preferentially occurring 
in regions of large-scale background temperature highs in a mantle with substantial 
viscosity increase with depth, although other models are possible. The effect of a 
relatively low viscosity asthenosphere is to have plumes neck down and attenuate 
there. The major mass anomalies causing the geoid highs associated with plumes 
appear to be in the lower mantle. 

Introduction: 

Linear seamount and island chains, such as the Hawaiian islands, have fre- 
quently been attributed to the passage of the lithosphere over deep convective upwel- 
lings (Wilson, 1963; Morgan, 1972; 1981). The age progression from the active 
“hotspot” to the guyots on the inactive end of the chain is particularly well 
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established for Hawaii (Jarrard and Clague, 1977; Dalrymple and Clague, 1976), and 
relative motion among the more prominent of these “hotspots” is constrained to be 
about an order of magnitude less than typical plate rates (Morgan, 1972; 1981; Enge- 
bretson et al., 1984; Chase, 1984). Therefore, the thermal plumes, or whatever process 
is responsible for hotspots, must be essentially stationary with respect to tectonic 
plate motions. 

The hypothesis of mantle plumes has not received universal acceptance, partly 
because much mid-plate volcanic activity is not easily associated with hotspot traces. 
For example, the Tertiary volcanic activity in eastern Australia (Pilger, 1982) and the 
recent volcanism near Easter Island (Bonatti et al., 1977) are actually “hot lines” 
rather than hot spot tracks. The Line Islands require either widespread contem- 
poraneous volcanism or several hot spot tracks (Schlanger et al., 1984; Epp, 1984b). 
Alternative explanations for mid-plate volcanism have usually involved propagating 
cracks or faults in the lithosphere (Betz and Hess, 1942; Turcotte and Oxburgh, 1973; 
1976; Sleep, 1974; 1984a; Solomon and Sleep, 1974), even though there is no resem- 
blance between the surface morphologies of mid-oceanic swells and other tensional 
features in the lithosphere such as mid-ocean ridges and continental rifts. 

The crack theory and the plume theory predict very different sub-lithospheric 
structures beneath a hot spot. These differences can be inferred by considering mid- 
plate swells such as the one associated with the Hawaiian Islands. These ~1000 km 
wide features are attributed to heating of the lower lithosphere as it passes over the 
hotspot (Detrick and Crough, 1978; Crough, 1978; Von Herzen et al., 1982; Epp, 
1984a). The topographic uplift appears to form within a few million years at the 
hotspot and then subside in a manner similar to young seafloor. The thermal origin 
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(within the lithosphere) of the swells is further indicated by their elevated heat flow 
(Von Herzen et al., 1982) and the systematics of volcano heights (Epp, 1984a). The 
formation of the hotspot swells is sufficiently rapid that bulk replacement of the lower 
lithosphere, as opposed to thermal conduction, is required (Detrick and Crough, 
1978). The replacement process could be intrusion of hot plume material into the 
lower lithosphere or, in the crack theory, bulk stoping or delamination of the lower 
lithosphere which then sinks as dense blobs into the underlying mantle. The two 
hypotheses therefore predict opposite types of structures deep beneath hotspots: a 
hot, low density plume or cold, sinking lithospheric material. 

Cases intermediate between these end members are conceivable, and it is neces- 
sary to clarify our terminology. By “plumes” we mean more or less cylindrical zones 
of upwelling with radii of the order of 10-100 km. Plumes might either be strong and 
supply the bulk of the heat needed to thin the lithosphere, or they might be weaker 
and act mainly as a trigger for delamination. Broad zones of mantle upwelling are dis- 
tinguished from narrow plumes. Mostly passive “blobs” in the upper mantle (Allegre 
et al., 1984; Batiza et al., 1984) which may cause chemical and isotopic anomalies in 
off-axis volcanism are also distinct from active plumes. We use the term “delamina- 
tion” to describe either thermally or mechanically triggered sinking of blobs of high 
viscosity material at the base of the lithosphere, i.e., convective instability. We distin- 
guish this process from lithospheric thinning due only to thermal erosion of the litho- 
sphere by a plume. 

Geophysical methods that might discriminate among these alternatives include 
modelling the gravity signatures of hotspot traces and studies of the deep seismic 
velocity structure beneath active hotspots. Seismic evidence would seem to favor the 



- 5- 


V 


plume hypothesis since certain hotspots such as Yellowstone are underlain by slow 
velocity material to a considerable depth below the lithosphere (Iyer, 1975; Hadley et 
al., 1976). Short-wavelength (<1000 km) gravity anomalies, although conspicuous, 
are largely the result of lithospheric thinning and compensated swell topography 
(Detrick and Crough, 1978; McNutt, 1984) and do not offer much direct information 
concerning dynamic processes deep in the mantle. However, very long-wavelength 
geoid anomalies (harmonic degrees <10) are relatively insensitive to contamination 
from lithospheric heterogeneity (Hager, 1983) and are most sensitive to the deep- 
seated density contrasts in the mantle (Richards and Hager, 1984, henceforth referred 
to as “RH”) that are presumably the result of convection. 

The general association of hotspots with long-wavelength geoid highs both glo- 
bally (Crough and Jurdy, 1980; Chase, 1979) and more locally (Kaula, 1970; Morgan, 
1972) suggests that, if hot plumes cause hotspots, topographic compensation dom- 
inates their geoid signature; otherwise, low density material would result in geoid 
lows. Paradoxically, subducted slabs representing cold, sinking material in the upper 
mantle are also associated with geoid highs (Kaula, 1970; Chase, 1979; Crough and 
Jurdy, 1980; Hager, 1984), indicating that dynamic surface topography is not the 
overwhelming effect there. Kaula (1970) explained this apparent contradiction as the 
result of differing rheology under hotspots and subduction zones. The effect of sub- 
ducting slabs encountering a relatively high viscosity lower mantle would be to reduce 
the resulting surface topography (“trench”) at long-wavelengths so that the positive 
geoid anomaly due to the dense slab dominates the geoid signature (Hager, 1984). 
Thus a conceivable alternative explanation for the geoid highs over hotspots is that 
cold, dense lithospheric blobs are delaminated at hotspots, resulting in geoid highs 
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just as do subducting slabs. 

The purpose of this study is to provide a test of the various hotspot theories by 
modelling the long-wavelength geoid anomalies with which they are associated. We 
use models of hot plumes and cold downwellings (slabs), which include temperature 
dependent rheology, to predict long-wavelength geoid anomalies which can then be 
compared with observations. We also consider other quantities of interest such as 
heat flow and long-wavelength dynamic topography, but they are more difficult to 
constrain via observation. A key question which we state now and expand upon later 
is: Could narrow mantle plumes be directly responsible for the geoid highs over 
hotspots, or are hotspots and, perhaps, plumes associated with more broad-scale tem- 
perature anomalies in the mantle which cause geoid highs? Emphasis in modelling 
will be placed on the conspicuous geoid high over Hawaii since it is the classic 
hotspot. However, in order to provide a more general observational base, we begin by 
analyzing quantitatively the relationship between hotspots and the geoid in a global 
sense, including information from recent seismic studies of mantle heterogeneity. 

Global Observations and Hotspots 

The Earth’s long-wavelength geoid (shown in Figure la referred to the hydros- 
tatic figure) has been well determined from observations of satellite orbits (Kaula, 
1963; Lerch et al., 1983), but its interpretation has remained somewhat enigmatic 
because of its lack of resemblance to surface features such as continents and mid- 
ocean ridges. However, by filtering out the lowest harmonics (degrees 2-3) which dom- 
inate the geoid spectrum (see Figure 2), it is obvious that many of the “intermediate” 
wavelength geoid highs are located oyer active subduction zones (Figure lb). This 
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subduction signal can be removed from the geoid with moderate confidence because of 
the high degree of formal correlation between slabs and the geoid at harmonic degrees 
4-9 (see Figure 3). Hager (1984) has presented a model, which we review in more 
detail below, that allows us to subtract a subduction geoid signal, calculated using a 
fluid dynamical model, from the observed geoid to obtain the residual geoid shown in 
Figure lc. 

The residual geoid is dominated by two large highs centered over Africa - north 
Atlantic and over the west-central Pacific. Crough and Jurdy (1980) and Chase(l979) 
recognized that residual geoid highs (left after subtracting slab effects) cover areas 
that include most of the world’s hotspots (marked with dots in Figure 1). Although 
not as striking as with all degrees included, upon filtering the lowest degree (2-3) com- 
ponents (Figure Id) we still find more “local” residual geoid highs over many of the 
hotspot provinces including Hawaii, Tasmania, Raton - Yellowstone - Bowie - Juan de 
Fuca, Christmas Island - Kerguelen - Crozet - Vema, Afar, Easter - Juan Fernandez, 
and Iceland - Madeira - Canary - Azores - Cape Verde - Rio Grande - Fernando. The 
Hawaiian anomaly is very striking and less likely than, e.g., Iceland to be contam- 
inated by plate boundary effects. There are conspicuous exceptions including Mt. 
Erebus, Samoa, and Bermuda which occur in pronounced residual geoid lows. Also, 
geoid highs remain over the Iranian-Caucasus-Tibetan highlands, which are related to 
convergence and thickening of the continental crust (Hager, 1983). 

We have selected our list of 47 hotspots (Table I) based on the compilations of 
Morgan (1981) and Crough and Jurdy (1980). Although exception may be taken with 
any of several inclusions or deletions (conceivably, only some hotspots are associated 
with plumes), this list probably represents the distribution fairly well. The 
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compilation by Burke and Wilson (1976) of 115 possible hotspots also exhibits a 
strong association with the low-degree geoid as shown by Crough and Jurdy (1980). 

The dynamic geoid response of the Earth to internal density contrasts depends 
quite strongly on the wavelength considered (RH). It is convenient, as well as instruc- 
tive, to calculate models for comparison to observations in the spectral domain using 
spherical harmonics. The spherical harmonic representation of hotspots which we use 
for statistical correlations is obtained by mathematically representing hotspots as 
point sources of equal (and arbitrary) strength on the surface of the Earth. We have 
made no attempt to selectively weight certain hotspots such as Hawaii, Iceland, and 
Kerguelen, which are surely more important than others, such as Raton, whose legi- 
timacy as hotspots may be questioned. 

The hotspot distribution spectrum is shown in Figure 2 along with the geoid and 
slab residual geoid spectra. These spectral amplitude plots are obtained from the 
square root of the sum of squares of harmonic coefficients at each harmonic degree /: 

A{l)= «/ £(c,l+s,l)/( 21 + 1 ) ( 1 ) 

V m =0 

The C{ m and are the cosine and sine coefficients for a fully normalized spherical 
harmonic expansion. The factor l/(2/-fl) is included because a random distribution 
of delta functions (hotspots) on a sphere will have a flat (“white”) spectrum with this 
normalization. The hotspot spectrum (Figure 2) is much whiter than either of. the 
geoid spectra; it is mildly peaked at degrees 1-2 with a striking peak also at degree 6. 
(Because the geoid is referred to the center of mass coordinate system, it has no 
degree 1 component.) 
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The correlation coefficient, r/, between the geoid and hotspots may be obtained 

from 


"b®/m ^lm ) 



where (c lm , are the geoid coefficients and (g/ m , /i /m ) are the hotspot coefficients. 
Cumulative correlations with several or many harmonic degrees simultaneously can be 
misleading since spectral power is not uniform (Eckhardt, 1984), so we examine only 
degree-by-degree correlations. Hotspots are significantly correlated with the observed 
geoid only at degree 2, but the low-degree correlations for hotspots vs. the slab resi- 
dual geoid are higher as shown in Figure 3. Confidence limit contours determined by 
a Student’s t test with 21 degrees of freedom are also shown in Figure 3. A 
confidence level of 0.95 implies that there is a 5% probability that the two sets of 
functions are random. 

The residual geoid is significantly correlated with hotspots at degrees 2, 4, and 6 
with some correlation at degree 3. Higher harmonics are essentially uncorrelated (no 
correlations are significant with >90% confidence for /= 7-20). The correlations at 
degrees 2 and 6 are particularly noteworthy because they correspond to peaks in the 
hotspot spectrum. Crough and Jurdy (1980) found a correlation coefficient of 0.85 at 
degree 2, significant with >95% confidence, that is even higher than our value of 
0.75; the difference arises from different methods of estimating the slab effects. 

The degree 2 correlation is made even more compelling by recent observations of 
seismic velocity heterogeneity in the lower mantle. Both the tomographic inversions of 
Clayton and Comer (1983) and the least-squares inversion of Dziewonski (1984) of P- 
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wave travel times show that low velocity in the lower mantle is very strongly corre- 
lated with low-degree (2-3) geoid highs (Hager et al, 1985). We also find that slow 
velocity and presumably hot, low density anomalies are well correlated with the 
hotspot distribution at degree 2 (r 2 =0.85). Figure 4 emphasizes this point by com- 
paring harmonic degree 2 maps of the slab residual geoid, a depth average of seismic 
heterogeneity from Clayton and Comer (1983), and the hotspot distribution. All of 
these fields closely resemble the entire low-degree residual geoid (Figure lc) because 
the geoid spectrum is so strongly peaked at degree 2. (The vertically averaged lower 
mantle P-wave velocity model is also peaked at degree 2 as shown in Figure 2.) These 
three phenomena are apparently related, and even though statistical correlations con- 
tain no information concerning cause-and-effect relationships, we form the following 
hypotheses: 1) The largest residual geoid highs are the result of long-wavelength topo- 
graphic highs that are dynamically supported by either broad scale or plumelike ther- 
mal anomalies; and 2) occurrences of hotspots (mantle plumes) are directly related to 
the broad-scale temperature structure of the lower mantle. That is, mantle plumes 
penetrating to the surface to form hotspots are preferentially located in regions of 
hotter than average mantle. The second hypothesis is suggested in part because the 
velocity anomaly spectrum inferred from tomography is redder than the hotspot spec- 
trum. 

Additional evidence comes from studies of upper mantle heterogeneity from sur- 
face wave studies (Masters et al, 1982; Woodhouse and Dziewonski, 1984; Nataf, 
Nakanishi, and Anderson, 1984; Tanimoto, 1986). At degree 2, there is a high velo- 
city feature in some models of the transition zone that correlates well with the geoid 
(e.g., Masters et al., 1982; Woodhouse and Dziewonski, 1984; Nataf et al., 1986). This 
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feature is even better correlated with subducted slabs than with the observed geoid 
(Hager, 1984; Richards and Hager, 1986) and appears to be associated with cold 
downwellings rather than hot plumes. Both the Woodhouse and Dziewonski and the 
Tanimoto studies show a remarkable correlation with both the residual geoid and 
hotspots at degree 6. Table II gives the degree 6 correlation coefficients, and the nega- 
tive signs indicate that slow shear velocity is correlated with both geoid highs and 
hotspots. Since there is a spectral peak in the hotspot distribution at degree 6, the 
correlation at degree 6 is expected if hotspots are to show a strong relationship to 
either the geoid or shear velocity anomalies. If the upper mantle is near the solidus, 
then shear waves could be very sensitive to elevated temperature. In Figure 5 we 
compare the degree 6 surface wave velocity heterogeneity, residual geoid, and hotspot 
distributions to illustrate the strength of an ~0.7 correlation coefficient at degree 6. 
(We should note that the lower mantle heterogeneity models do not correlate 
significantly with either hotspots or the geoid at harmonic degrees >4. Lack of reso- 
lution may be at fault. Also note that degree 6 is the one harmonic degree for which 
slabs do not correlate well with the observed geoid. However, it is only when the 
dynamically modelled degree 6 slab geoid is removed that the hotspot distribution 
shows good correlation with the geoid.) Figure 5 shows that the degree 6 hotspot peak 
represents the large groupings of hotspots (e.g., Christmas Island - Kerguelen - Crozet 
- Vema) rather than individual spacings which appear to be random. (The hotspot 
spectrum beyond degree 10 is essentially white.) Again, we can formulate a testable 
hypothesis concerning hotspots: Plumes are directly related to either heating or chem- 
ical heterogeneity in the upper mantle at degree 6, which may be a dominant 
wavelength for their formation. 
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These observations suggest that we formulate a quantitative global test to deter- 
mine whether mantle plumes might be directly responsible for the density contrasts 
that cause the large-scale residual geoid features as well as the seismic velocity 

anomalies. The question then becomes that of whether reasonable models of 

hotspots, either mantle plumes or delaminating lithosphere, can explain the slab resi- 
dual geoid. The alternative, of course, is that hotspots are only symptoms of a 
broad-scale thermal field or, perhaps, compositional heterogeneity in the mantle. 

The correlations of Figure 3 suggest to us that much of the long-wavelength slab 
residual geoid is causally related to hotspots. We assume that for each harmonic 
degree this relationship can be written in the linear form 

residual geoid — ( dynamic response)* (hotspot distribution) + (noise) 
or, for example, 

( c /m ) = b l ( 9lm Am ) + ( ™ise ) (3) 

(Note that this assumption is in contradiction of our second hypothesis, which 

assumes an additional component of heterogeneity, the “long-wavelength back- 
ground.”) From our analysis we obtain the least-squares estimates for the dynamic 
response functions, , shown in Figure 6. The coefficients (<7/ m A m ) are in units of 
hotspots, and the spectral “response” curve is in the rather peculiar units of 
geoid/hotspot. Although we cannot reliably determine the response at degree 5, 
where the correlation is poor, it appears that the response is a relatively smooth, 
monotonic function of harmonic degree, consistent with a dynamical filtering process 
as discussed below. Also shown are the values of b t obtained if, instead of assuming 
that all the error in estimation is the result of other density heterogeneity signals in 
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the geoid, we perform the mutual correlation of residual geoid and hotspots under the 
pessimistic assumption of equal noise in each signal. The extra noise on the geoid 
(left) side of equation 3 may be due primarily to mismodelling of the subduction sig- 
nal; at these wavelengths the geoid can otherwise be considered to be perfectly meas- 
ured. The “equal noise” response is not substantially different than the initial model 
(equation 3), so our response curve is robust at least in this respect. Unfortunately, 
the least-squares fits for the coupling coefficients, 6/ , are less well constrained as 
shown by the 2<r error bars at each harmonic degree. The best fitting response ampli- 
tudes give about a factor of 8 decrease from harmonic degree 2 to degree 6. This spec- 
tral shape is largely that of the residual geoid, since the hotspot spectrum does not 
show the same long-wavelength bias. 

In addition to the global association of hotspots with geoid highs, we can also 
use the local ~13 m geoid high (degrees 4-10) centered on Hawaii to constrain our 
models (see Figures lb,c). The contours are not elongated in the direction of the older 
islands and seamounts in the chain (toward the northwest), so it is difficult to explain 
this long-wavelength signal as an effect of the lithospheric swell itself. Intraplate den- 
sity contrasts do not contribute substantially at these wavelengths (Hager, 1983). 
The anomaly actually appears to be elongated toward the upstream direction 
(southeast), suggesting, perhaps, that the active Hawaiian shield is lagging slightly 
behind a deep thermal source. 

The spectral content of the Hawaiian geoid anomaly is difficult to assess quanti- 
tatively because it is necessary to arbitrarily select some spatial subdomain within 
which to perform spectral analysis. However, Figures 7a, b show that about 10 m of 
the 13 m signal occur in the harmonic degree 4-6 band, while less than ~3 m occur in 
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the degree 7-12 band. Figure 7c shows that the degree 10-20 geoid signal over Hawaii 
is almost zero, and also verifies the lack of any consistent correspondence between 
geoid highs and hotspots (noted above) in this wavelength band. (Note, however, 
that there is a strong shorter-wavelength signal over Yellowstone.) Since the Hawaiian 
swell is of relatively small width (~1000 km), the lack of degree 7-20 signal makes it 
an implausible source for the longer-wavelength positive geoid anomaly. These obser- 
vations for the isolated case of Hawaii are consistent with the pronounced long- 
wavelength bias of the inferred global hotspot geoid response curve of Figure 6. We 
use both the local and global observations to discriminate among long-wavelength 
geoid responses for the competing hotspot models discussed below. 

Dynamic Response Functions 

At this point we review some basic ideas about how long-wavelength geoid 
anomalies are generated in a viscous, convecting planet like the Earth. Chase and 
McNutt (1982) and Hager (1983) have shown that only about 20 meters out of a total 
long-wavelength geoid signal of about 200 meters can be generated by compensated 
topography and lithospheric or crustal thickness variations, e.g., the geoid high over 
the Tibetan Plateau. Therefore, most of the geoid must result from the internal den- 
sity contrasts that drive convective flow: subducted slabs, mantle plumes, or broader 
scalelength variations. 

Interior density contrasts drive flow that causes deformations of the surface, the 
core-mantle boundary, and possibly, internal chemical boundaries. At very long 
wavelengths (/ <10) the lithosphere has effectively no long-term flexural strength 
(McKenzie and Bowin, 1976), and deformation will occur rapidly compared to the 
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timescale for convection (RH). These deformed surfaces have an important effect on 
geoid anomalies. In order to correctly model the long-wavelength geoid, a fluid 
dynamic Earth model must be used to calculate the geoid contributions due to these 
boundary deformations. It has been shown by many authors (e.g., Pekeris, 1935; 
Runcorn, 1964; Morgan, 1965; McKenzie, 1977; Parsons and Daly, 1983; Ricard et al, 
1984; RH) that dynamic compensation due to boundary deformation is of dominant 
importance in determining the geoid. Since induced boundary deformations cause 
geoid anomalies that are of opposite sign and comparable magnitude to the geoid due 
to interior density contrasts, long-wavelength geoid anomalies are the difference of 
large numbers. The details of boundary deformation depend strongly on the viscosity 
structure of the mantle, so the geoid is a sensitive indicator of mantle structure (RH). 

If the viscosity structure varies only radially (i.e., is spherically symmetric), then 
a given density contrast Sp lm {r) at radius r excites only an lm th harmonic flow field 
and causes only lm lh harmonic boundary deformation. Since solutions for linear 
(Newtonian), spherically symmetric viscosity may be superposed, we can obtain the 
total harmonic geopotential coefficients from 

U lm = ^p^-jGl( r ) S Plm{ r ) dr (4) 

where 7 is the gravitational constant, R the Earth’s radius, c the core radius, and 
G t (r ) is the dynamic response function or kernel. This kernel is independent of the 
azimuthal order m and contains contributions from both boundary deformations and 
the density contrast itself. In RH we showed how to analytically calculate G)(r ) for 
spherically symmetric, incompressible, self-gravitating Earth models. 


& 
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Response functions for both whole mantle flow and chemically layered flow are 
shown in Figure 8 with lower/upper mantle viscosity ratios of 1, 10, and 100. Note 
that although flow velocities depend on the absolute value of viscosity, the stresses, 
boundary deformations, and geoid depend only on the relative values. Free-slip boun- 
dary conditions are imposed at the core and at the surface; the difference between 
no-slip and free-slip is discussed in RH. For uniform viscosity and whole mantle flow 
(Figure 8a) the geoid response is always negative because of the overwhelming gravi- 
tational effect of the deformed upper boundary. Decreasing the viscosity of the upper 
mantle causes less deformation of the upper boundary (Figure 8b, c) and tends to 
drive the geoid response toward more positive values. Therefore, both the size and 
magnitude of the geoid response are strongly affected by relatively mild changes in 
viscosity with depth. These pressure induced changes can occur either gradually due 
to compaction or abruptly due to phase changes; phase changes probably do not oth- 
erwise strongly affect the flow field (Richter and McKenzie, 1981). However, a chemi- 
cal discontinuity acting as a barrier to radial flow will deform and also affect the 
geoid. This forces the geoid response to zero at the boundary (in the same way that 
we get perfect compensation at the surface and core) and generally reduces the mag- 
nitude of the response functions (Figure 8d,e,f). 

In addition to the response function G)(r) we have also calculated dynamic 
impedance functions, z*(r) which give the ratio of induced topography to observed 
geoid (RH). Unfortunately, the Earth’s dynamically supported topography is obscured 
at long-wavelength due to gravitationally compensated continental masses and ther- 
mal plate thicknesses, and there is no consistent correlation between topography and 
gravity. At the present time we cannot reliably estimate the global long-wavelength 
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(degrees<10) dynamic topography associated with mid-oceanic swells or hotspots. 

An obvious application of the response functions of Figure 8 is in modelling the 
subducted slab geoid signal (see Figure lb). A dynamical model using deep seismicity 
to locate subducted slabs was developed by Hager (1984). By associating slabs with 
approximately 0.1 gm/cm 3 density contrast and convolving these mass anomalies 
with various response functions, the following conclusions were reached based on com- 
parison of the observed and predicted geoids: (1) The magnitude and sign of the 
observed geoid response is consistent with a relatively small degree of dynamic com- 
pensation at the surface with most of the dynamic compensation occurring at the 
core-mantle boundary. (2) The sign of the response function is positive for harmonic 
degrees 2-9 in the upper mantle, requiring a viscosity increase with depth. (3) Chemi- 
cally layered models require about a factor of 5 more density contrast associated with 
subducting slabs than expected. (4) The best-fitting two-layer viscosity model is that 
of whole mantle flow with a factor of about 30-100 increase in viscosity between the 
upper and lower mantle (see Figure 8c). 

Figures 8b,e,f show that it is possible to have positive responses in the upper 
mantle and negative ones in the lower mantle for a variety of models. This at first 
suggests one easy solution to the problem of geoid highs over both slabs and plumes, 
namely, that the mass anomalies associated with plumes are primarily lower mantle 
features. (More realistic models that include a low viscosity asthenosphere and high 
viscosity lithosphere added to model U10 give more positive upper mantle kernels but 
maintain negative kernels in the lower mantle.) This idea turns out to be basically 
correct in the numerical plume models which are discussed below. However, we feel 
that it is important to assess the impact on our response kernels caused by neglecting 
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the large viscosity variations expected to be associated with slabs and plumes. One 
motivation for numerical modelling is the hypothesis that these viscosity variations 
might be responsible for the paradox of having geoid highs over both slabs and 
plumes. 

An alternative explanation is that the geoid highs over hotspots are due to 
delamination of cold lithosphere, consistent with the slab results. With the upper 
mantle “calibrated” by the geoid response of cold, subducting slabs, it is straightfor- 
ward to estimate the geoid signature of unstable lithosphere sinking below a hotspot 
if we can estimate the amount of high density material present. We present both 
types of models in the following sections. 

Delaminated Blobs 

Although the uppermost 30 km of the lithosphere under Hawaii behaves elasti- 
cally (e.g., Watts, 1978), the portion of the thermal lithosphere below the elastic layer 
should behave as a cold, dense, high viscosity boundary layer. It is convectively 
unstable and might sink into the mantle, or “delaminate.” For Hawaii, if all of the 
swell topography is attributed to delamination, the flux of delaminated blobs is com- 
parable to that of slabs at subduction zones. The age of the crust around Hawaii is 
about 90 Myr and the apparent thermal age after the ~ 1 km uplift of the swell is 
around 25 Myr (Epp, 1984a). The elevation and hence the average mass anomaly in 
the lithosphere is proportional to the square root of age, so the delamination is 
equivalent to subduction of a 20 Myr plate at the rate of hot spot migration (100 
mm/yr) across the ~1000 km width of the swell. We assume that large-scale hor- 
izontal motions in the mantle are much less than plate velocities (e.g., Hager and 
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O’Connell, 1979), so that the locus of delamination is roughly fixed with respect to 
the mantle, not the plates. If we assume 10 Myr as a characteristic time of transit 
through the upper mantle, then, for a slab of material 1000 km long (along swell) by 
1000 km wide with an excess mass per unit area of 3.3z 10 6 kg/m 2 (associated with 
the uplift of the swell), the excess mass of lithospheric material in the mantle beneath 
Hawaii is ~3z 10 18 kg. 

This load can be convolved with the preferred geoid response curves (“U100”) 
for subducted slabs (Figure 8c) to estimate the long-wavelength geoid. If the excess 
mass is roughly distributed in a cylinder <1000 km in radius and 1000 km deep 
beneath Hawaii, we predict about an 6.3 m geoid high over Hawaii for harmonic 
degrees 4-9, which is about half the observed signal. However, the spectral response 
for this model is not nearly as strongly peaked at the lowest degrees as the global 
response curve (see Figure 18 for a comparison) owing to the small horizontal scale of 
the load. The degree 2-10 topographic downwarp is about 85 m, a value not likely to 
be resolved by analysis of bathymetric data given other perturbing influences. The 
lithospheric swell itself will generate relatively little geoid signal since it is isostatically 
compensated at shallow depth; whatever signal is generated will also be essentially 
’’white” at low harmonic degrees since the swell is only ~1000 km wide. 

The shorter wavelength (1 > 10) geoid contribution could be as much as 10 m if 
there is little compensation, but this value depends strongly on the detailed viscosity 
structure of the upper mantle (RH). It is evident that the short- wave length geoid is 
not elevated 10 m near the hotspot, and about 300-400 m of downwarp is required to 
keep the / >10 geoid anomaly small. If this downwarp occurs, a place on the seafloor 
would be expected to experience subsidence (superimposed on the lithospheric swell 
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due to delamination) as it approached the hotspot, and then uplift as it drifted 
further west. Since we cannot predict the detailed timing of the delamination it is 
difficult to model the upstream (east) side of the hotspot. On the west side of Hawaii, 
simple thermal contraction should produce about 350 m of down warp in the first 11 
Myr if the lithosphere is reset to a thermal age of 25 Myr. If this downwarp is super- 
imposed on an uplift of ~300-400 m due to rebound as the lithosphere moves away 
from the sinking, delaminated blobs, we obtain approximately a neutral net 
uplift/subsidence on the downstream side of the hotspot. Although this computation 
is crude, this might explain the fact that the Hawaiian swell has not substantially 
subsided as far as 10° West of the hotspot as shown in Figure 9. (This topographic 
variation could conceivably be attributed to greater activity of the hotspot 11 Myr 
ago.) This description of swell topography due to delamination can be contrasted with 
the following model based on heat flux from a mantle plume. 

Plume Kinematics 

Before describing our numerical plume models, we consider a simplified 
kinematic description of a plume for Hawaii that relates the heat flow, mass flux, 
temperature excess, and swell topography. Assuming for now that there is no trig- 
gered lithospheric delamination, the plume must diverge widely enough beneath the 
lithosphere to account for the swell width and must also supply enough heat to thin 
the lithosphere. That is, the flux of positive buoyancy from the plume should equal 
the rate of production of positive buoyancy in the swell. Using, as before, a migration 
rate of 100 mm/yr, a width of 1000 km, an elevation of 1 km and a density of 3300 
kg/m 3 for the swell, the net flux m of negative buoyancy is 10 Mg/sec. The actual 
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mass flux in the plume is the buoyancy flux divided by the fractional density contrast 
in the plume, mp/6p. For thermal expansion the density contrast is 6p=paAT , 
where AT is the excess plume temperature and a is the volume thermal expansion 
coefficient, 3z 10 -5 /°K. The resulting heat flux from the Hawaii hotspot, Cm /a, 
(C = l.2x 10 3 J/kg °K is the specific heat) is 4.2x 10 11 W or about 1% of the global 
mantle heat flux. The volume of flow, Q, through the plume necessary to make the 
swell is inversely proportional to the temperature contrast 

Q = m /apAT (5) 

(For reference, a plume with Ar=100°K, 100 km diameter, and m=10 Mg/sec 
ascending as a cylindrical plug will have an ascent velocity of 4 m/yr.) 

Next, the flow from the plume must diverge widely enough to produce the 
observed swell. To a first approximation this flow can be considered to be the super- 
position of radial flow from the plume through an asthenospheric channel and the 
horizontal drag produced by the motion of the plate over the hot spot. The vertically 
averaged velocity in the asthenosphere far away from the center of the plume is 

y plume =(Q/27rrA)a r (6) 

where r is the horizontal distance from the plume, A is the thickness of the astheno- 
sphere, and a r is a unit radial vector from the hotspot. The vertically averaged velo- 
city in the asthenosphere from drag at the base of the plate is about half the plate 
velocity if the base of the asthenosphere is moving much more slowly than the overly- 
ing plate: 

v 0 = {Q /2irrA )k r -v L /2& x (7) 


where the x direction is positive upstream from the plume and v L is the plate 
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velocity. The upstream stagnation point occurs at r s = Q/ttAv l . The stagnation 
streamline assuming a 100 km asthenosphere thickness is shown for temperature con- 
trasts of 300°K and 1000°K (Figure 9) and the flux computed above. The 1000°K 
curve is narrower than the swell, but the 300°K curve is a fairly good fit to the edge 
of the swell, which extends about 500 km ahead of the hotspot. (The topographic 
contours for the Hawaiian swell in Figure 9 have been adapted from Schroeder, 1984, 
who computed the anomalous seafloor topography in the Pacific ocean after correc- 
tions were made for isostatic loading, sediment thickness, seafloor age, etc. Note that 
the islands and seamounts themselves involve many kilometers of topography, but 
these loads occur at short wavelength and are partially supported by lithospheric 
flexure, e. g., Watts, 1978.) 

The preceeding calculation is crude, but it shows that the swell shape and buoy- 
ancy flux are kinematically consistent with a plume delivering several or many hun- 
dreds of degrees excess temperature and also with a reasonable limit for the heat flux 
due to the Hawaii hotspot. (This heat flux is probably at least a factor of 2 or 3 too 
high for an average hotspot, since 47 hotspots would otherwise account for half of the 
Earth’s entire heat budget.) It remains to be seen whether a mantle plume fitting this 
surface kinematic model can also satisfy the geoid observations and constraints from 
convection theory on the thermal structure of plumes. We consider more refined 
models in the next section. 
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Fluid Dynamical Models of Plumes 

Our basic idea of a plume is that of a narrowly confined, stable, hot upwelling 
from a deep boundary layer due to bottom heating. The seismic D' ' layer at the 
bottom of the Earth’s mantle may be in part the result of heating of the base of the 
mantle by radioactive, specific or latent heat within the core. Most estimates for the 
amount of core heat flux are somewhat less than 10% of the total geothermal flux 
(Gubbins et al, 1979). If there is a plume under Hawaii it is probably less than ~200 
km in radius as evidenced by the width of the volcanic trace (Morgan, 1972b). This 
small dimension is consistent with a very high effective Rayleigh number and the 
strong temperature dependence of viscosity, as shown by the stability analysis of 
Yuen and Peltier (1980). 

The thermal structure of a plume is dependent upon such unknowns as the rheo- 
logical laws of the mantle and plume material, the amount of heat being vented, the 
superadiabatic temperature drop, whether the plume is chemically distinct from the 
upper mantle, and the possible influence of partial melting. Thus, we do not know 
the plume structure very well. By contrast, we can make an educated guess as to the 
thermal structure of a subducted slab or even a delaminated lithospheric blob. Our 
approach is to address some general questions about plume dynamics and compare 
the results to the observations and hypotheses discussed above. In particular, we ask: 
(1) How does the low viscosity of the plume alter the geoid responses derived for 
purely layered viscosity? (2) Is the thermal buoyancy of a narrow plume consistent 
with both the observed amplitude of hotspot geoid anomalies and reasonable limits on 
the amount of heat flow due to a hotspot? 
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The main difficulty in modelling is the extreme horizontal variations in viscosity 
expected for thermal plumes. No analytical methods exist to treat this problem in 
detail, so a numerical solution is required. We need to calculate very long-wavelength 
(~10,000 km) stress fields to obtain long-wavelength surface deformation and geoid 
estimates, but we must also resolve the strong short-wavelength (• — 10-100 km) plume 
structures which drive the flow. The finite element method is well suited to this prob- 
lem because of the advantage of variable grid size, and a typical grid for our plume 
models is shown in Figure 10. In all of the calculations that follow we have used 10 
km horizontal grid spacing for treating the prescribed buoyancy forces and viscosity 
variations for both plumes and slabs, thus allowing good resolution of thermal struc- 
tures as thin as 50 km or less. Much smaller spacings are possible by further packing 
the element lattice at the origin, but this was not necessary except as a check on solu- 
tion accuracy. We have also used 20 km vertical spacing in the top 200 km and bot- 
tom 100 km of the mantle to properly resolve viscosity changes in the 
lithosphere/asthenosphere and core-mantle boundary (D ' ' ) region. Rotational or 
reflection symmetry imposed about r =0 (r is the radial distance from the plume 
center) or x =0 gives a total effective width of 10,000 km in both cylindrical and 
Cartesian geometry. (In the cylindrical calculations we assume a free-slip boundary 
at the outer boundary of the cylinder, r=d.) The numerical code is based on a 
penalty method formulation (Hughes et al., 1979) of the Stokes flow problem (steady, 
incompressible, very low Reynold’s number flow with spatially variable Newtonian 
rheology). The code handles nonlinear rheology by damped iteration upon the viscos- 
ity field. 
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We are currently limited to numerical solutions in two dimensions (2-D Carte- 
sian or cylindrical geometry with axial symmetry) rather than the spherical geometry 
for which we earlier showed analytical solutions. The two-dimensional results that fol- 
low are represented in the horizontal spatial wavelength domain, just as we 
represented our spherical, analytical models (Figure 7) in spherical harmonics. In 
Cartesian coordinates (appropriate for subducted slabs) we have the approximate spa- 
tial wavelength equivalence, 

\ F « 2 *R /Vl{l+l) 

where \p is the Fourier transformed spatial wavelength, / is the corresponding spher- 
ical harmonic degree and R is the Earth’s radius. In cylindrical coordinates (r ,z) 
with no 6 dependence, we use the Fourier- Bessel transform of the spatial coordinate r 
(see Appendix A). Therefore, for axial symmetry (appropriate for plumes) we have the 
approximate wavelength equivalence, 

2nd/p n ' 2nR /Vl(l + 1) 

where p n 1 is the n th zero of the derivative of the zeroth order Bessel function, 
J 0 (r ), and d is the radius of the cylindrical domain. For our numerical grid, d =5000 
km and =10,000 km for Cartesian geometry, so the maximum allowable 
wavelengths in both cases correspond approximately to spherical harmonic degree 4. 
(Due to memory limitations, this was the maximum horizontal dimension our com- 
puter model could handle with accuracy and still provide high resolution at the center 
of symmetry.) The first and second “overtones” in both coordinate systems 
correspond to harmonic degrees 8 and 12, respectively. 
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We now consider three different cases in which the low viscosity of a plume may 
affect the geoid signal. The first case is that of hot blobs (such as those investigated 
experimentally by Olson and Singer, 1985) guided in their ascent by a narrow, pre- 
established pathway. The upper surface deformation and geoid due to these solitary 
blobs can be adequately modelled by the analytical theory (RH). Unless the blobs are 
very closely spaced, i.e., connected, there will not be an effective low viscosity stress 
pathway to the upper surface (or the core-mantle boundary), so induced surface 
deformations will not be very different from those calculated using the ambient man- 
tle viscosity structure. Their buoyancy will act approximately as point sources, with 
respect to the long-wavelength flow-stress field, embedded in a high viscosity back- 
ground (mantle). The particular case of a hot blob impinging on the lithosphere is a 
special case which we will discuss later. A more closely spaced string of low viscosity 
hot blobs might behave more like the plume models described next. 

The other two plume types with which we are concerned are like the more classi- 
cal, steady-state structures in which vertical flow is very rapid with respect to mantle 
flow as a whole and is nearly uniform, with convergence at the bottom toward the 
plume and divergence at the top. The probable narrowness (<200 km) of mantle 
plumes implies that the zones of flow convergence/divergence will be of equally small 
dimension (perhaps that of the D ' ' layer and the asthenosphere); we temporarily 
ignore these complications. We can envision two end-member classes of these steady- 
state plumes for which the low plume viscosity may affect the geoid signature. Mantle 
plumes may be essentially low viscosity “pipes” in which low viscosity rising fluid is 
contained by the relatively rigid walls (mantle). The other type is more typical of 
mantle convection models in which the radial excess temperature profile decays 
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smoothly away from the axis of the plume (e.g., Parmentier et al., 1975; Yuen and 
Schubert, 1976; Boss and Sacks, 1985). If the plume is not distinct chemically from 
the surrounding mantle, a long-lived plume should evolve to this latter state as the 

"s 

surrounding mantle is heated conductively and, possibly, by viscous dissipation. We 
concentrate on this “thermal halo” case in the models that follow, returning to the 
“pipe” case when a comparison is needed. 

The thermal profile for the halo model can be parameterized by a characteristic 
width, r 0 , and by a peak (axial) excess temperature, AT 0 , which should be essentially 
independent of depth if the plume rises nearly adiabatically. The temperature profile 
near the plume axis must satisfy the condition 

4^ = 0 at r =0 (8) 

dr 

Since most of the flow and transport will occur very near the axis due to temperature 
softening of the mantle material, the dynamics of the plume is sensitive to the excess 
temperature profile. In accord with equation 8, in some models we have adopted the 
form used by Loper and Stacey (1983) 

AT = AT„(l-(r/r 0 ) 2 l (9) 

where A To is the peak excess temperature in the plume and r 0 is a characteristic 
plume radius. More commonly, we modify (9) to include a thermal “halo” surround- 
ing the mobile near axis region by a temperature profile of the form 

AT = AT 0 exp[-(r /r 0 ) 2 ] (10) 

where r 0 gives the characteristic thermal width of the plume. Estimated values for 
AT 0 associated with the superadiabatic increase across D ' ' vary widely up to a 
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maximum of about 1000°K (Verhoogen, 1973). This parameter is not well constrained 
by observation, and we consider a wide range of values in our models. 

The density contrast in the plume is given by 

8p = -p 0 aAT (11) 

where p 0 is the density at the background mantle temperature and a is the volume 
coefficient of thermal expansion. We have taken p 0 = 5.14 g/cm 3 and Q'=3a: 10" 5 /°K. 



depth, so a temperature contrast of only 100°K results in an excessive uplift of 4.3 
km. Therefore, viscous drag must limit flow in the plume. 

The viscosity of mantle minerals is a strongly decreasing function of tempera- 
ture, and we have used the exponential form 

V = Vo™P (-0AT/T 1 ) (12) 

which also closely mimics the form used by Loper and Stacey (1983) for a characteris- 
tic temperature T' =2300°K and /?<35. For example, with AT =800°K and 0=35, 
we obtain more than five orders of magnitude viscosity decrease from the colder sur- 
rounding mantle to the hot plume axis. Stress-dependent rheology, by diffusing 
viscous stresses away from the plume, lowers the effective value of 0 (Christensen, 
1984). 

Note that for large values of /? or AT 0 the combination of exponential depen- 
dences in expressions 10 and 12 will restrict most of the flow to a very narrow region 
near the axis, even though most of the thermal buoyancy, or “halo,” may lie outside 
of this region. Expressions 10 and 12 allow us to characterize a wide variety of 
dynamical behavior by using 0, A T 0 , and r 0 to specify, independently, the maximum 
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viscosity contrast, the thermal buoyancy, and the characteristic width of the plume. 
This parameterization could closely match most published temperature and viscosity 
profiles for plume models. 

We calculate geoid response models for a series of plume geometries of increasing 
complexity illustrated in Figure 11. (Again, in these models, we assume a thermal 
structure, rather than obtaining a completely dynamically consistent one using a con- 
vection code.) In our first model (A) we consider a simple plume of characteristic 
width r 0 = 70 km and peak excess temperature 700°K through the entire depth of the 
mantle. Using expressions 10-12 to specify the load and laterally varying plume 
viscosity structure, we have calculated the total geoid response by obtaining the 
induced boundary deformation from a finite element solution. Strictly speaking, the 
depthwise, wavelength dependent response functions of equation 4 do not exist for 
laterally varying viscosity; the different spatial wavelengths are mutually coupled. 
However, we can construct “pseudo response functions” for comparison with analytic 
kernels by considering all of the geoid at a particular wavelength to be due only to 
the load component at that wavelength. The depthwise pseudo response functions 
for a plume with about 4 1/2 orders of magnitude viscosity contrast (/?=35) are 
shown in Figure 12. 

In the absence of lateral viscosity variations, the flow and stress fields due to a 
density contrast of a given spatial wavelength are independent of disturbances at 
other wavelengths. Two-dimensional geoid response functions similar to in equa- 
tion 4 can be calculated analytically as a function of spatial wavelength, and the dot- 
ted lines in Figure 12 show responses for a uniform mantle at spatial wavelengths 
corresponding to harmonic degrees 4, 8, and 12. These functions are exactly the same 
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in Cartesian and cylindrical geometry for a given wavelength (see Appendix A). The 
response is uniformly negative due to the dominance of boundary deformation, and 
the curves are similar to those for the Ul model of Figure 8. 

The curves for the low viscosity plume are more negative in the upper mantle 
than for the uniform viscosity case in the upper mantle, but the effect is rather small 
in comparison to that of depthwise viscosity variation or chemical stratification (Fig- 
ure 8). One reason for this is that much of the plume buoyancy, or thermal halo, lies 
outside of the very low viscosity plume core and the resulting change in surface defor- 
mation is relatively small. The change that does occur causes the geoid response to 
become more negative due to enhanced long-wavelength boundary deformation. We 
note here that changes in the outer radius boundary condition on the cylindrical 
domain (at r =d ) from free-slip to rigid have less than a h% effect upon the geoid 
calculations. Therefore, the finite cylinder radius probably has less effect than other 
neglected effects such as the dynamical interaction of plumes. The long-wavelength 
deformation field is even less sensitive, so it is safe to conclude that the 5000 km max- 
imum radius for the finite element grid is not a severe limitation in these calculations. 
Extrapolation to a degree 2 wavelength may be more questionable, but in that case 
the effects of self-gravitation (RH), sphericity, and self-compression (Ricard et al, 
1984; Hong and Yuen, 1985) are even more important. The response curves in Figure 
12 are truncated at 200 km depth, because the solution accuracy degrades (>l% 
error in the l =4 surface deformation compared to analytical solutions) for loads 
above this level. 

The narrowness of a single plume gives essentially equal loading at all 
wavelengths of interest (a spectrally “white” load) so that the relative geoid response 
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at each wavelength can be obtained by integrating along the response curves. From 
Figure 12 we see that the low plume viscosity causes enhancement of the shorter 
wavelength (/ =8,12) geoid and surface deformation compared to the longer 
wavelengths (/=4). The low plume viscosity allows more efficient transfer of the 
buoyancy forces (“head”) in the plume to the upper and lower surfaces, so that the 
load is effectively closer to these boundaries. The net result is that the low viscosity 
associated with plumes tends to make the geoid response spectrum more blue. 

The geoid responses in Figure 12 are for a particular choice of parameters. To 
determine the sensitivity to parameter values, we conducted a range of numerical 
experiments. The magnitude of the effects of ft upon the total geoid signature is 
shown in Figure 13 for a wide range of the parameters ft and AT 0 . (The relative size 
of this effect also increases with the plume radius as demonstrated for the pipe models 
in Figure 16.) Our numerical experiments show that for (ftA T 0 /T' )<6 (less that 2 
1/2 orders of magnitude viscosity contrast) the geoid signal is enhanced by a rela- 
tively modest factor of 20% or less. Unless the temperature dependence of viscosity is 
much stronger than given by the rather high value of ft=35, a plume temperature 
excess of at least 600-800°K will be required to substantially affect the geoid signature 
of this type of mantle plume. 

The plume radius, r 0 = 70 km, in the example of model A was chosen so that the 
amplitude of the long-wavelength geoid signal (/ =4,8) would be comparable to that 
observed over Hawaii (~13 m) for ATo = 700 0 K- The geoid elevation (per °K tem- 
perature contrast) from this halo model with no viscosity perturbation {/3=0) is 
(0.72, 0.72, 0.55) m/100°K for / =( 4, 8, 12) scale lengths. These values can be scaled 
by the cross-sectional area of the plume (~To 2 ) within a few percent accuracy over 
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the range 30 km < r 0 < 100 km. The Ageoid curves of Figure 13 give the geoid sig- 
nal (/?>0) in excess of the unperturbed signal (0=0). The size of the geoid anomaly 
will scale roughly as AT 0 r 0 2 unless the effects of low viscosity are very large. 

In matching the observed geoid signal, we must not exceed the approximate 
v upper limit of ~10 Mg/sec buoyancy flux (see previous section). Figure 14 shows the 
buoyancy flux, normalized to a background mantle viscosity »7 0 =10 21 Pa-sec, for 
r 0 =70 km plumes with varying Ar o and (3. For /3=35, AFo^^OO 0 ^ we obtain 
~200 Mg/sec buoyancy flux, so the mantle viscosity must be raised by a factor of 20 
to 2x 10 22 Pa-sec to lower the flux to 10 Mg/sec. (Flow velocities in all of our examples 
scale inversely as j/ 0 .) Flow in these plume models is limited by viscous drag in the 
surrounding mantle; the buoyancy flux varies little with depth except very near the 
top and bottom of the mantle, and the values in Figure 14 are calculated at mid- 
depth. Figure 15 shows that the buoyancy flux scales linearly with the 4 lh power of 
the radius, r 0 4 , as expected for flow in a long, narrow conduit. 

We now consider a “pipe” model in which the temperature profile is steplike, 
i.e., 

A T = constant , r <r 0 (13) 

A T = 0, r >r 0 

Both temperature and viscosity are uniform within the pipe. Figure 16 shows the 
geoid effect of low viscosity for a suite of pipe models parameterized by the radius, r 0 , 
and by the viscosity contrast relative to the background mantle viscosity. The unper- 
turbed geoid signal (pipe viscosity = 1.0) is (0.63, 0.63, 0.48) m/lOO^K for ro=70 km 
(almost the same as for the halo model) and scales as AT r 0 2 as for the halo model. 
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We get about 10 m of geoid signal for 800°K excess temperature in the pipe. 

Now we ask: How does the pipe model compare to the halo model in perturbing 
the geoid signal? We begin by locating a pipe model, r 0 =70 km, in Figure 16 that 
gives about the same perturbation (~30%) to the 1—4 geoid signal (see Figure 13) as 
did halo model A. This requires a viscosity contrast of about a factor of 0.006 within 
the pipe. The buoyancy flux for the pipe models is mapped in Figure 17 and scales as 
the square of the excess temperature (buoyancy x force). For a temperature contrast 
of 800°K, radius r 0 =70 km, and a pipe viscosity contrast of 0.006, we obtain a buoy- 
ancy flux of 36 Mg/sec, or about a factor of 3 1/2 more than for the comparable halo 
model. Therefore, given the restriction on the maximum allowable buoyancy flux and 
given the requirement for the amplitude of the geoid signal, the pipe model is not as 
efficient as the halo model in perturbing (increasing) the size of the geoid signal from 
the uniform viscosity value. Similar comparisons show that this difference between 
the pipe and halo models persists, to a varying degree, for other plume radii. Also 
shown for comparison in Figure 17 (dashed lines) are the fluxes calculated analytically 
for an infinitely long pipe in a rigid mantle, but with the same pipe viscosity and 
buoyancy forces; i. e., classical Poiselle flow (proportional to the fourth power of the 
pipe radius). For very low pipe viscosities (.001-.0001 mantle background), the flow in 
numerical experiments is limited by the finite pipe length. For high pipe viscosity (.1- 
1.0 mantle background), induced flow in the mantle becomes significant and the flux 
in numerical experiments is higher than the analytical result. 

In the models that follow we find that large perturbations in the geoid signal 
(e.g., sufficient alone to account for geoid highs over both subducted slabs and 
hotspots) often require unacceptably large buoyancy fluxes. Therefore, in seeking to 
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understand the largest effects on the geoid of low plume viscosity, we concentrate on 
the halo model. As we discussed previously, many of our inferences about mantle 
viscosity structure are based on our analysis of geoid anomalies associated with sub- 
ducting slabs (Hager, 1984; Hager and Richards, 1984). Having seen the effects of low 
viscosity plumes on geoid response curves, it behooves us to address the related prob- 
lem of how the high viscosity associated with subducted slabs affects these responses. 
Before presenting more sophisticated models for hot plumes, we investigate the geoid 
signature for a very simple model of a cold subducting slab. 

We model a high viscosity slab by considering (in Cartesian geometry) a 100 km 
wide slab, density contrast 0.1 g/cm 3 , and having two orders of magnitude higher 
viscosity than the surrounding mantle. The “slab” for this first example extends 
throughout the depth of the mantle, so this case and the previous plume models 
might be taken to simulate the rising and falling hot and cold columns in a very high 
Rayleigh number, bottom heated, whole mantle convection ■ system. The pseudo 
response curves for the slab model in Figure 12 show that for /=4 and 8 as well as 
for / =12 in the upper mantle, the high slab viscosity causes the response to be much 
less negative (more positive) than for a uniform mantle viscosity. If the background 
viscosity layering in the mantle gives a “marginal” upper mantle response, perhaps 
with a zero crossing such as for model U10 (Figure 8b), both slabs and plumes resid- 
ing entirely in the upper mantle are qualitatively consistent with positive geoid 
anomalies. 

We now repair some of the obvious inadequacies in plume model A. In model B 
we have added a high viscosity lithosphere, a low viscosity asthenosphere or outlet 
channel, and a low viscosity D' ' layer which simulates the lower boundary layer 
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feeding the plume (see Figure lib). The plume parameters are similar to those of 
model A with Ar 0 =700°K, /?=35, and a slightly diminished radius, r 0 =60 km. The 
total resulting buoyancy flux, 11.0 Mg/sec, accounts for the maximum heat flux for 
Hawaii if the mantle reference viscosity, q 0 , is 10 22 Pa-sec. The peak velocity in the 
plume is then 2.3 m/yr. 

The long-wavelength dynamic topography (/ =4,8) for plume model B is about 
165 m. The total long-wavelength geoid signature (1=4,8) is 10.9 m, and the degree 
4 and 8 responses are plotted along with the estimated global response curve in Fig- 
ure 18. (The /= 4 and / =8 geoid anomalies from the plume calculation are divided 
by 4 as plotted in Figure 18; in a real Earth the power concentrated at these 
wavelengths due to the limited radius of the finite element geometry will actually be 
smeared over an interval of approximately 4 harmonic degrees.) Although we have 
obtained a geoid amplitude and buoyancy flux compatible with observation, there are 
three serious problems with this plume model which are very instructive: 

(1) The background mantle viscosity of 10 22 Pa-sec required to keep the buoy- 
ancy flux down to a reasonable value is probably too high a viscosity to assign to the 
entire mantle based on post-glacial rebound and rotational data (Peltier, 1981; Yuen 
et al., 1982). (This viscosity is probably not too high for the lower mantle, and we 
shall explore this point further in model C.) 

(2) The geoid spectrum from model B is much too white (Figure 8). If we sum up 
contributions for all wavelengths less than 500 km, we get a total geoid signature of 
40 m and 2.6 km of dynamic uplift (in addition to the uplift due to lithospheric thin- 
ning). Such a model is clearly not acceptable for any hotspot, including Hawaii. 
Increasing or decreasing the viscosity of the lithosphere in the model does not 
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substantially alter the large dynamic uplift; it is mainly the result of allowing the 
outlet of a strong narrow plume to be within 200 km of the surface. 

(3) Subducted slabs in the upper mantle part of model B will not give a positive 
geoid signal, even if we assign the same viscosity to a downgoing slab as to the litho- 
sphere (see Appendix B). In order to get positive upper mantle geoid response curves 
similar to those in model U100 of Figure 8 (required in order to fit the observed geoid 
anomalies over subduction zones), the viscosity of the asthenospheric channel must be 
at least a factor of 10,000 smaller than the underlying mantle, which we find implau- 
sible. 

Problem (l) in model B can be eliminated by decreasing the average viscosity 
contrast in the plume by simply lowering /? or by simultaneously lowering A T Q while 
increasing r 0 (in order to conserve the total geoid signal). Problems (2) and (3) are 
more difficult and are addressed in model C. 

We have yet to find a model that can explain the apparent geoid signatures of 
both slabs and plumes. An obvious way to approach this problem is to start with a 
reference model whose upper mantle responses are neither strongly negative nor posi- 
tive. Model C (Figure 11) satisfies this criterion by including a high viscosity litho- 
sphere, a low viscosity zone extending to 400 km depth, a moderate viscosity increase 
through the transition zone, and a higher viscosity lower mantle. We have also 
included a low viscosity D' 1 layer in which the viscosity profile is determined by a 
temperature gradient of 7°K/km (/?=35) in accord with the boundary layer model of 
Stacey and Loper (1983). This layer does not significantly affect the lower geoid har- 
monics or mid-mantle buoyancy flux, but it is included for consistency with the idea 
of plumes originating at the core-mantle boundary. If there are no lateral viscosity 
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variations, the response functions for this model (Figure 19, dotted lines) are small 
and negative in the lower mantle and are small and generally positive in the upper 
mantle. 

To obtain the geoid response for model C we have used the same plume parame- 
ters as in models A and B for the lower mantle except for the plume radius, r 0 , which 
we have increased to 100 km to compensate for the generally smaller low-degree 
response functions. In the two low viscosity upper mantle layers the plume must 
diminish in radius if the plume buoyancy flux is approximately constant throughout 
the entire mantle (i.e., steady-state). If we assume that the temperature excess at the 
plume’s center is also constant (both the plume and mantle are adiabatic) and that /? 
does not vary with depth, then the plume radius, r 0 , should decrease as the fourth 
root of the layer viscosity as demonstrated above. A constant flux plume of 100 km 
radius in the lower mantle should neck down to ~32 km radius in the upper mantle 
since r)—>rj 0 / 100. Since the plume buoyancy and geoid signature are proportional to 
r 0 2 , the upper mantle plume will contribute relatively little to the total long- 
wavelength geoid signature. Most of the geoid signal will result from the lower mantle 
plume. The reduction of width as the plume rises will also be enhanced if a large 
fraction of the plume partially melts and causes a large decrease in the plume viscos- 
ity. For a temperature excess of 300-700°K the plume might encounter the solidus as 
little as ~1000 km above the core-mantle boundary (Anderson, 1981). Plumes in the 
upper mantle may only be streamers of melt from more substantia] solid state plumes 
in the lower mantle. 

The depthwise pseudo response functions for this plume model are shown in Fig- 
ure 19. The lighter line in the upper mantle emphasizes that even though the 
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normalized response is of the same order for the upper and lower mantle, the upper 
mantle buoyancy multiplying this response (equation 4) is very small by comparison. 
Again, the response curves are more negative than for the purely layered case, but 
when integrated through the lower mantle the resulting geoid spectrum is much 
stronger at the longer wavelengths (/ =4) than at shorter wavelengths (/ =12). Note 
that this would be the case regardless of whether or not the low viscosity in the 
plume affects the response (compare the solid and dotted curves of Figure 19). The 
plume radius of 100 km for this mode! results in a geoid signature of similar ampli- 
tude to the global response curve as shown in Figure 18. The geoid spectrum is now 
much more “red” (biased toward long-wavelengths) than for our other plume models 
and is a more satisfactory approximation to the global response spectral shape than 
any of our previous models. (The total signature for l— 4,8 is less than 5 m, so scaling 
the plume radius up to about r 0 =140 km will account for the Hawaiian anomaly.) 
The short wavelength geoid signature (l >12) is only a few meters, so model C does 
not suffer the problems of very large geoid anomalies at shorter wavelengths that we 
found for model B. The shear stress, r rz , at the base of the lithosphere drops from 
about 10 bars at a distance of 200 km from the plume center to only about 1 bar at a 
distance of 1200 km, so the plume will not drive much plate motion. The shear stress 
at plume center is over 200 bars and might result in erosion of the lithosphere. 

The buoyancy flux for model C (r 0 =100 km) is 85 Mg/sec for a lower mantle 
reference viscosity i? 0 =10 22 Pa-sec, however, so unless the average viscosity of the 
lower mantle is as high as 10 23 Pa-sec, this model must be considered unacceptable. 
But if we use a weaker temperature dependence for effective viscosity (likely due to 
stress-dependent effects) or a lower excess plume temperature, we can greatly reduce 
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the flow while preserving the geoid signature. For example, a reduced value of (3=22 
produces essentially the same geoid signature as for /3=35 with a buoyancy flux of 
only 10.5 Mg/sec for rj 0 =10 22 Pa-sec. Alternatively, for AT 0 =300°K and 35, 
essentially the same geoid response may be produced by a plume of radius 200 km in 
the lower mantle and with a buoyancy flux of only 0.52 Mg/sec. A wide range of 
plume parameters can, therefore, match the geoid and flux constraints either globally 
or for the particular case of Hawaii. 

Both the low viscosity plume geoid response curves and the purely layered 
viscosity curves (dotted) in Figure 19 will give the correct long- wavelength spectral 
characteristics for hotspots. The overall negative response functions (resulting in geoid 
highs for low density plumes) in the lower mantle overwhelm the upper mantle plume 
signal because the plume’s radius is smaller in the low viscosity upper mantle. This 
necking down effect may be very important to the dynamics of mantle plumes, espe- 
cially plumes of lower mantle origin, and appears to be a neglected phenomenon in 
both numerical and experimental modelling. 

To emphasize the relative importance of vertical stratification of mantle viscos- 
ity, we have included two additional models in the spectral response plot of Figure 18. 
First, model C ' (/?=0) is the same as model C, except that there is no viscosity con- 
trast between the plume and surrounding mantle. Its spectrum is somewhat more 
“red” than for /?= 35, and, according to the numerical experiments (model A type) 
discussed above, this purely layered model is probably sufficiently accurate for up to 
two orders of magnitude viscosity contrast within the plume. We can also perform 
this type of calculation ((3= 0) analytically for spherical Earth models including all 
harmonic degrees. For a 300°K plume of 200 km radius (no viscosity contrast) in the 
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model C type layered mantle, we obtain the response given by dashed lines in Figure 
18. (The slight differences between this curve and model C ' illustrate the effects of 
sphericity.) This response curve gives a reasonable, although not optimal, fit, at least 
for / <8, to the observed global response curve. Clearly, increasing mantle viscosity 
with depth can result in a strong low-degree bias in the geoid signature of a plume (or 
any convective upwelling or downwelling), regardless of the viscosity of the plume 
itself. This is one of the main points to be gained from our numerical models. 

Our geoid models involve the balance between dynamic surface deformation and 
the (plume) load. Figure 20 shows the actual radial profiles of surface deformation 
from models B and C. Shown for comparison are Hawaiian swell profiles A-A ' and 
B-B ' from Figure 9, which cross the island of Hawaii and a point 500 km “down- 
stream,” respectively. Figure 20 shows the relative richness in long-wavelength sur- 
face deformation in model C compared to model B. Dashed lines also give the defor- 
mation with no viscosity contrast in the plumes (/?— 0). For our “preferred” model C, 
the predicted long-wavelength deformation outside of the swell is several hundred 
meters, which is probably not resolvable from bathymetric anomalies. However, the 
swell topography itself is not explained by this model. Furthermore, since the ridge- 
like swell topography extends far toward the WNW from the active hotspot, lithos- 
pheric thinning (effectively resetting the thermal age of the lithosphere) is obviously a 
more satisfactory explanation (Detrick and Crough, 1978). The predicted short- 
wavelength (<500 km wide) deformation from model C of more than 1 km will be 
reduced by lithospheric flexure and masked to a great extent by the ~8 km seamount 
topography (volcanic edifice) itself, which does not appear in Figure 20. 
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We now have an acceptable working model for a plume derived long-wavelength 
geoid. The main virtues of model C are that it has a generally negative lower mantle 
geoid response for l <8 and that it involves a substantially reduced upper mantle 
contribution. We must now determine if this model is compatible with the observed 
subducted slab geoid response. 

Contrary to what we found for model A, for model C, with a 100 km wide slab 
with the same viscosity as the lithosphere and extending to 1100 km depth, the 
response (Figure 19, heavy-dashed line) is now more negative than for the purely lay- 
ered model and produces almost a null geoid signature for slabs at degrees 4 and 8. 
The slab load is coupled relatively more efficiently to the lithosphere than before, 
resulting in more surface deformation and a more negative geoid response. 

This coupling may not occur in the real Earth since the lithosphere at subduc- 
tion zones is weaker than normal (e.g., Sleep, 1979; Hager and O’Connell, 1981). If 
we simulate the weak plate boundary by reducing the effective viscosity of the litho- 
sphere by two orders of magnitude within 100 km of the subducting slab, the slab is 
supported more by the high viscosity of the lower mantle. This results in less long- 
wavelength surface deformation and a much more positive response (Figure 19, heavy, 
dash-dot line). 

Comparing the slab-with-weakened-lithosphere response with either the plume 
response or the purely layered response, we see that it is possible to have geoid highs 
over both plumes and slabs. This in itself is not too surprising considering the results 
of model A. However, the detailed rheology of subduction zones is of even more 
importance than the influence of low plume viscosity. The “calibration” of the upper 
mantle geoid response using subducted slabs depends strongly on the 
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parameterization of the lithosphere. Hager’s (1984) conclusions, based on a simple 
two-layer parameterization of rheology, that the slab-geoid observations favor pene- 
tration of the 670 km discontinuity by slabs and require a one to two order of magni- 
tude increase in mantle viscosity through the upper mantle are valid, but the details 
of the upper mantle structure are not well resolved. In Appendix B we have included 
a more systematic analysis of models for rheological variations associated with sub- 
ducting slabs. On the basis of these models we conclude that the lithospheric weaken- 
ing effect is the most efficient way to make the slab geoid response more positive than 
for purely layered viscosity. A more extensive treatment of the large-scale dynamic 
support of subducting slabs, including such effects as dip angle and non-linear rheol- 
ogy, might yield important new insights, but such a study is beyond the scope of this 
paper. 

In order to obtain average subducted slab velocities of about 100 mm/yr in 
model C, it is necessary to set the lower mantle reference viscosity t] q to 10 22 Pa-sec. 
This value, which is substantially higher than recent estimates from post-glacial 
rebound (Peltier, 1981) and rotational data (Yuen and Sabadini, 1982), served as the 
reference viscosity in our model C plumes. We can think of two resolutions of this 
apparent paradox. One is that the rotational response is actually on the high viscosity 
branch (O’Connell, 1971). The other is that postglacial rebound samples transient 
rheology while convection responds to steady-state rheology. 

Finally, we consider an upper mantle plume for model C, i.e., one that originates 
above 670 km depth. In order to get the 1 — 4,8 purely layered responses to become 
negative in the upper mantle, it is necessary to have more than 3 orders of magnitude 
viscosity contrast in the plume. As before, a plume radius >70 km is required to 
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produce the observed geoid anomaly over Hawaii. However, we again obtain too much 
heat flow (even for an upper mantle viscosity as high as 10 21 Pa-sec) just as we did 
for plume model B. This problem results from the low viscosity of the plume, and 
this is the main difference between our calculation and Morgan’s (1972) estimate of a 
75 km radius upper mantle plume. The plume can only overcome the effects of 
viscosity stratification (required by slabs) by stronger coupling of buoyancy from the 
lower parts of the plume to the surface through the low viscosity channel. It is very 
difficult then to produce large positive geoid anomalies over primarily upper mantle 
plumes without violating reasonable limits on plume heat flow. Since the conditions 
under which model C can satisfy the slab observations are somewhat extreme (high 
viscosity slab and very weak lithosphere), this conclusion is difficult to avoid by con- 
trivances in the upper mantle viscosity structure. This is not to say that strictly 
upper mantle plumes do not exist, but only that they are an unlikely source of long- 
wavelength geoid highs. 

Larger-scale Plumes in the Upper Mantle? 

The final issue we address is that of more broad-scale heating associated with 
hotspots in the upper mantle, or, perhaps, a solitary blob of hot material impinging 
on the base of the lithosphere. (The correlation of surface wave velocities with the 
degree 6 geoid and hotspot expansion may be symptomatic of a hotter than average 
asthenosphere.) In order to fit the observed long-wavelength bias described above, the 
basic horizontal scalelength of heating surrounding a hotspot must be >5000 km. 
This will lower the average viscosity of the upper mantle on a scale much larger than 
the depth of the upper mantle, resulting in geoid response curves that are more 
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positive , not more negative. In other words, coupling of the buoyancy forces to the 
lithosphere will be weakened by higher temperature and lower than average viscosity, 
resulting in less surface deformation. Therefore, given the requirement of increasing 
upper mantle viscosity with depth (based on slab modelling), we cannot explain geoid 
highs over such broadscale hot blobs in the uppermost mantle. Again, this does not 
exclude the existence of hotter than average asthenosphere associated with hotspots 
(which we would expect at the head of any thermal plume), but our models imply a 
deeper source for the associated geoid anomalies. 

In this same vein, it is curious that the strong correlations of the residual geoid 
with hotspots at degrees 2 and 4 are not accompanied by significant correlations 
between hotspots and upper mantle seismic velocity variations (Richards and Hager, 
1986). This also suggests a deep source associated with, if not caused by, plumes, and 
is consistent with the greater sensitivity of the lower geoid harmonics to density 
heterogeneities at great depth in the mantle. Lastly, we note that the surficial evi- 
dence for mantle plumes indicates that they are very narrow, at least in the upper 
mantle. The Hawaiian swell is only of order 1000 km in width, and the volcanic 
shields (islands) form a much narrower track within the swell. An active thermal 
source at the base of the lithosphere an order of magnitude wider seems unlikely. 

Discussion 

The geoid models we have presented show that there is no inherent contradic- 
tion in having geoid highs produced by both subducted slabs in the upper mantle and 
mantle-wide plumes. If the thermal buoyancy of plumes directly generates the 
observed geoid highs, we can offer some restrictions on their characteristics: 
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(1) Plumes that are of uniform strength in terms of thermal buoyancy (approxi- 
mately constant in radius), or that are primarily upper mantle phenomena, cannot 
produce the very long-wavelength bias (degrees 2-6) of the observed geoid. They are 
also unlikely candidates to produce the more local (/ >4) geoid highs over hotspots. 

(2) Plumes that are of radius ~100-200 km (A T 0 Pa200-700°K) in the lower 
mantle, but that become much more confined as they rise due to decreasing mantle 
viscosity, can produce the low-order geoid (degrees 2-6). However, such plumes must 
be restricted to about 3 orders of magnitude or less viscosity contrast in order not to 
exceed maximum estimates for heat flux from the core. Much narrower, low viscosity 
plumes (e.g., Loper and Stacey, 1983) can also account for the heat flux, but they 
cannot cause the long-wavelength geoid anomalies. 

(3) More local, intermediate wavelength (/ >4) geoid anomalies associated with 
hotspots may be caused by either the plume itself or triggered lithospheric delamina- 
tion due to the plume. Unfortunately, this study has not revealed a method based on 
geoid models to discriminate between these alternatives. It is remarkable that the 
delamination could indeed cause large (>5 m) geoid anomalies, but this process is not 
able to produce the low-order (/ =2-6) observed geoid without producing much more 
pronounced and consistent local (/ >6) anomalies over hotspots. Of particular impor- 
tance is the temperature of the sublithospheric upper mantle, which would be heated 
by plumes but cooled by delaminated blobs. (The mantle at lithospheric depths is 
heated by both processes.) The strong association of hotspots with both low velocity 
upper mantle (degree 6) and lower mantle (degree 2) as revealed by seismic studies, 
along with the relatively stationary nature of hotspots with respect to plate motions, 
implies that they are at least symptomatic of deep thermal processes. None of these 
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observations supports passive lithospheric delamination (e.g., due to cracks) as the 
primary source of either hotspots or the observed geoid. 

The correlation of hotspots with low-order residual geoid highs (Figure lc) and 
with seismic velocity anomalies (Figure 4) does not necessarily imply that the thermal 
anomalies within the plumes associated with hotspots are the primary cause of these 
phenomena. The considerable effort we have made to explore the conditions for which 
this is possible should not be taken to imply that we necessarily believe in such a 
strong role for plumes. Even our “preferred” model C appears unsatisfactory in some 
respects. In order to explain the low-degree (2-3) geoid, nearly all 47 of our selected 
hotspots would have to be associated with very large plumes (e. g., radius=200 km 
and average temperature contrast 300°K), and with conspicous “local” (/ >4) geoid 
anomalies. Hawaii and a few other major hotspots might fullfill this prescription, but 
many others would seem unlikely candidates. 

Further insight into this question may be gained from consideration of the global 
hotspot distribution. Hotspots are distributed almost randomly over about half the 
Earth’s surface (contained by the large residual geoid highs), but are almost absent in 
the other half of the globe. Stefanic and Jurdy (1984) have claimed less than 1 % 
likelihood that such a broad-scale bimodal pattern is random, and we suspect that 
the location of hotspots is controlled by some other global thermal pattern in the 
mantle that limits their surface access to the areas of the large-scale geoid highs. An 
alternative explanation related to a process of delamination is not evident. 

A candidate control mechanism is found by considering the reconstruction of 
paleo-subduction zones proposed by Chase and Sprowl (1983). If, as they claim, the 
Pangean supercontinental assemblage was surrounded by subduction zones at ~125 
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Myr B.P., then the major residual geoid lows (Figure lc) correspond to areas which 
have experienced subduction as the American continents have swept westward during 
the last ~100 Myr. These are also areas largely devoid of hotspots, with some excep- 
tions such as Yellowstone. Intense shearing in the mantle and/or thermal quenching 
due to deeply subducted material may block or completely shut off mantle plumes. 
Strong plumes such as Hawaii, Iceland, and Kerguelen have probably been shielded 
from subduction for 100 Myr or more. Chase and Sprowl also point out that other 
hotspots such as Mt. Erebus, having only recently escaped a subducting slab, are rela- 
tively weak newcomers. 

This conceptual model provides an alternative explanation for the correlation of 
low seismic velocity at very long wavelength in the lower mantle with both the largest 
geoid highs and hotspots. High seismic velocity in the lower mantle may be giving us 
a broad-scale map of dead slabs that have been deposited in the lower mantle during 
the last 100-200 Myr, and these positive mass anomalies, combined with negative 
lower mantle response kernels, would generate geoid lows. Both Figures 8 and 19 
show that this would not be a contradiction of the fact that currently subducting 
slabs in the upper mantle cause geoid highs. Lower mantle response functions are 
probably negative even with a low viscosity asthenosphere. We can roughly estimate 
the degree 2 geoid from ancient subducted slabs in the lower mantle as follows: 

Subducted slabs currently residing in the upper mantle cause about +20 m of 
degree 2 geoid for an average upper mantle response of +0.3 (see Figure 8c). For a 
purely layered such as used for model C 1 model we calculate an average lower man- 
tle degree 2 response (using the spherical, analytic model of RH) of about -0.1. If we 
assume that the ancient subduction zones have put 10 times as much lithosphere into 
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the lower mantle over the past 100 Myr as currently resides in the upper mantle due 
to recent subduction, then we estimate ~70 m amplitude degree 2 geoid lows associ- 
ated with these old subduction zones. This is about the right size to explain the 
current low-degree geoid. Also, dead slabs in the lower mantle will produce a very 
attenuated signal at higher harmonic degrees, resulting in a geoid spectrum strongly 
peaked at degree 2. Further reddening of the geoid would result from shearing and 
diffusing away of short wavelength slab heterogeneities. Since at least 70% of the 
Earth’s heat flow is involved in the cooling of lithospheric plates (O’Connell and 
Hager, 1980), this explanation is satisfactory in that the largest geoid features are 
related to the dominant mode of convection. This scenario, with a relatively passive 
lower mantle, is similar to the recently proposed whole mantle convection model of 
Loper (1985), although we have been motivated by different observations. The man- 
tle convecting system, dominated by internal heating and driven by subduction of the 
upper boundary layer (plates), would be semi-transparent to mantle plumes resulting 
from a relatively small heat flux from the core. Such plumes could, of course, be 
much smaller in radius than those required to directly cause the long- wavelength slab 
residual geoid, which might result from long-wavelength variations in background 
temperature. 

If the Pangean supercontinent was stable for a long period of time, the geoid 
may have a long-term memory of that episode as proposed by Anderson (1982). A 
ring of subduction around this supercontinent may also have resulted in antipodal 
rifting and broad-scale upwelling in the central Pacific which is also “remembered” by 
a large residual geoid high. Hotspots shielded from subduction beneath the super- 
continent would have caused long-term heating beneath the continental lithosphere. 
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This may have eventually led to the breakup of the stable configuration as evidenced 
by the many hotspots along the mid-Atlantic ridge and the African continent. Thus, 
the two convecting systems interacted strongly at this point. The chemically buoyant 
supercontinent may have stabilized a degree 2 convection mode on Earth of which the 
present geoid is a fossil, and the correlation of hotspots with the low-degree geoid is 
symptomatic rather than causal. 

This hypothesis stands in contrast to the possible dominant role of hotspots on 
Mars and Venus. The largest gravity and topography anomalies on Mars are due to 
the Tharsis bulge (Sleep and Phillips, 1985) which contains huge shield volcanoes. On 
Venus, the large gravity highs over the highland areas are most easily attributed to 
large mantle plumes in a relatively isoviscous mantle (Kiefer et al., 1986). A central 
question that remains, in our view, is that of just which density contrasts in the 
Earth’s mantle cause the large-scale shape of both the geoid and the lower mantle 
velocity variations. We cannot completely resolve this issue on the basis of our study, 
even though only a restricted class of hotspot models can account for the geoid. 
However, because not all hotspots are associated with local geoid highs, even our 
“preferred” plume model is not a very satisfactory explanation for low-degree (2-3) 
geoid highs. More broad-scale sources appear to offer a less problematic explanation. 

This work might be improved by substituting a variety of fully dynamic convec- 
tive plume solutions for our generalized plume models. We hope that some of our 
obvious oversimplifications may spur other workers to predict long-wavelength geoid 
anomalies from their plume models. Improvement in the resolution of seismic velocity 
anomalies in the deep mantle can be expected, since this is a relatively recent area of 


research. 
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Appendix A: Dynamic Geoid Anomalies in Cylindrical Geometry 

The governing first-order differential equations for incompressible flow at very 
low Reynold’s number include the continuity condition 

V-v = 0 


the equations of motion 


V-r + 6pg = 0 

and the Newtonian constitutive law 


t = -p I 4- 2 rje 


where v is the velocity vector, g the gravitational acceleration, dp the density con- 
trast, r the stress tensor, e the strain rate tensor, jj the viscosity, p the pressure, and 
I the identity matrix. In cylindrical coordinates (r ,z) we have 
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Tgz — P+ 2 V dz 
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dv r 
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We can eliminate the r dependence by an appropriate Fourier-Bessel respresen- 
tation. For example, for a finite cylinder of radius a the radial velocity must vanish 
at a , so we can write the inverse discrete transform 

Mr,*)- E ■’,”(* P, I — ) = F f' [»") 

n =1 i a J 

where [t> r n ] is the discrete Fourier-Bessel transform of v T (r ,z ) and p n 1 is the n lh zero 
of the first-order Bessel function J With some foresight concerning the use of ortho- 
normal expansions, and with F 0 respresenting the zeroth order transform, we also 
write: 

1 K] 

Tzz = Fo 1 \t"} 

Trz = F^ \t”) 

P =Fo l [p n > 

6 p = F 0 - 1 \ 6 p n } 


The stresses t 69 and t tt will turn out to be mixtures of J 0 and J 1; but will be 
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formally eliminated from the equations in a later step. 

Some useful orthogonality and completeness properties of Fourier-Bessel 
transforms are summarized below (see Sneddon, 1951). 

Orthogonality: 
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where p m and p q can be either the zeros or the zeros of the first derivative of . 
Representation (completeness): 
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Useful derivatives: 
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for any real k . 
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By taking JrJ 0 (kr) [eg. (1)] dr and using the expansions for v z and v r , we 
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obtain a transformed equation: 
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where D =— and the wavenumber is k— . Similarly we obtain for Equations 6 
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and 3 
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If we use equations 1 , 
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from equation 

2, we obtain 
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Equations l' , 6' , 3' , and 2' are identical in form to the two-dimensional Cartesian 
formulation of Cathles (1975) and Hager and O’Connell (1981) 
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where we have dropped the wavenumber superscript n . This system of equations can 
be solved analytically for v and t via a propagator matrix technique for arbitrary 
layering of viscosity with depth ( z ) (see Hager and O’Connell, 1981, for examples). 


Thus the depth dependence for cylindrically symmetric flow in layered media 
can be solved in exactly the same way as for 2-D Cartesian coordinates. For a 


prescribed cylindrically symmetric loading problem driven by bp 


Pn 


we can write 
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with the Cartesian solution driven by 5p(k), where k=2ir/\ and X is the spatial 
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wavelength, if 



a 


In other words, cylindrical and Cartesian flow are “wavenumber equivalent”. Note 
also that in an effectively infinite domain (a — »-oo) we can use the integral transform 

00 

v r( r , z ) = J v r ( 2 , a ) J l(<* r W a 
0 

to again obtain equations 8. In practice, we have used the discrete transform to 
represent our finite element results in the wavenumber domain. 

Appendix B: Slab Geoid Models with Variable Viscosity 

We have tested two series of models. The first (“LlOO”) has a 100 km thick high 
viscosity (100%) lithosphere with a uniform background mantle viscosity (%). Geoid 
anomalies are calculated in a box of width 4000 km (corresponding to 1=5 with 
reflection symmetry imposed at the left side boundary). The finite element grid is the 
same one shown in Figure 10. The 100 km thick slab extends to 1100 km depth (Fig- 
ure Bl) and is assigned a density contrast of 0.1 g/cm 3 . (The density contrast actu- 
ally normalizes out of the response curves.) Pseudo response functions are calculated 
at four depth intervals ranging from 200 to 1100 km depth (accuracy is not sufficient 
for loads above 200 km) and are plotted in Figure B2 for the first and second har- 
monics (/ =5,10) of the box. 

Model LlOO.u (Figure B2, dotted lines) gives the purely layered response. In 
model LlOO. slab (solid lines) the slab is assigned viscosity 100%. In models 
Ll00.slab200 and Ll00.slab400 (long and short dashes) the lithospheric viscosity is 
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lowered to 7? 0 within 200 km and 400 km, respectively, of the center of the downgoing 
slab. As shown in Figure B2, the effects of these rheological variations are relatively 
mild. Stronger coupling of the slab load to the lithosphere through the high viscosity 
slab causes the responses to become more negative in the deeper part of the slab. 

The second series of models (“LU100”) is the same as the L100 series except that 
the background viscosity of the lower mantle (below 670 km) is increased to 100tj 0 . 
Figure B3 gives the pseudo response curves for models LUlOO.u, LUlOO.slab, 
LUl00.slab200, and LUl00.slab400 (dotted, solid, long-dashed, and short dashed lines, 
respectively). Here the effect of decreased lithospheric strength near the slab is more 
pronounced. The slab is less strongly coupled to the lithosphere and more strongly 
coupled to the high viscosity lower mantle. This results in less upper surface deforma- 
tion and gives a much more positive geoid signature for models LUl00.slab200 and 
LUl00.slab400. Again, we cannot resolve the geoid anomalies accurately for loads 
above 200 km depth. However, our parameterization is probably too coarse to 
represent the subduction process in this complicated zone even if a finer grid spacing 
were possible. 
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Table Captions 


Table 1 - Hotspot locations used to obtain the spherical harmonic expansion of 
the hotspot distribution. 


Table 2 - Correlation coefficients at harmonic degree 6 between upper mantle 
shear velocity models (Tanimoto, 1986; Woodhouse and Dziewonski, 1984) and 
the observed geoid, the slab residual geoid, and the hotspot distribution. 
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TABLE I: Hotspot Locations 


Lat . 

Long. 

Name 

50 

7 

E i f e 1 . Belg i um 

23 

6 

Hoggar Mountains. Algeria 

21 

17 

Tibe s t i , Ch a d 

13 

24 

Jebel Marra, Sudan / Darfur 

4 

9 

Mt. Cameroon 

- 3 

36 

Lake Victoria / East Africa 

12 

42 

Afar / Ethiopia 

- 1 2 

44 

Comores Islands 

-21 

56 

Reuni on 

-45 

45 

Cr o z e t 

-45 

65 

Kerguelen 

- 35 

80 

Christmas Island. Indian Ocean / Amsterdam 

-40 

150 

Tasmania 

45 

- 1 1 1 

Yel lovstone 

37 

-104 

Raton. New Mexico 

27 

-113 

Baja California / Guadalupe Seamount 

53 

- 135 

Bowie Seamount / Kodiak Seamounts 

46 

-128 

Juan de Fuca / Cobb Seamount 

20 

-155 

Hawaii 

- 29 

- 1 40 

MacDonald Seamount / Cook--Austral Islands 

- 18 

- 146 

Mehetia / Society Islands / Tahiti 

- 27 

- 1 29 

Pitcairn Island / Cambier Islands 

3 

167 

Caroline Islands 

- 1 1 

- 1 39 

Marquesas Islands 

- 27 

- 109 

Easter Island 

0 

- 91 

Galapagos Islands 

- 27 

-80 

San Felix, Nazca Plate 

- 34 

-79 

Juan Fernandez, Nazca Plate 

72 

-8 

Jan Mayen 

64 

- 20 

Iceland 

30 

-60 

Ber mud a 

38 

- 28 

Azores 

33 

- 17 

Madeira 

28 

- 17 

Canary Islands 

29 

- 29 

New England Seamounts / Great Meteor Seamount 

15 

-24 

Cape Verde 

-4 

- 32 

F ernando 

- 17 

-25 

Arnold Seamount 

-21 

-29 

Trindade 

-8 

-14 

Ascension 

-16 

-6 

St . Helena 

-37 

-12 

Tristan de Cunha 

-42 

0 

Discovery Seamount 

-54 

4 

Bouvet 

- 32 

16 

Verna Seamount 

-78 

167 

Mt . Erebus 

- 15 

-168 

Samoa 



-64- 


TABLE II: Degree 6 Correlations (significance) 

Upper Mantle Velocity Model 

TSVf350kml TSVf250km) W&Df200-500km) 

geoid -.63(>98%) -,67(>99%) -.58(>95%) 

residual geoid -.72(>99%) -.74(>99%) -.65(>99%) 

hotspots -.54(>95%1 -.49f>90%) -.70(>99%) 
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Figure Captions 


Figure 1 - (a),(b) Observed long-wavelength geoid (Lerch et al., 1983) referred to 
the hydrostatic figure of the Earth (Nakiboglu, 1982). In (b) the degree 2-3 
components have been removed. Hotspot locations are indicated by black 
dots. Continents and plate boundaries are also shown. Geoid lows are 
shaded; cylindrical equidistant projection. 

(c),(d) Residual geoid after the subducted slab geoid model of Hager 
(1984) is removed. 

Figure 2 - Log-log comparison of root mean square harmonic coefficient ampli- 
tudes as a function of spherical harmonic degree. Units are as follows: 
Observed and residual geopotentials, /R (fraction of geopotential at sur- 
face), (7 is the gravitational constant, M the mass of the Earth, and R the 
Earth’s radius.) Hotspot distribution, 4.1zl0 6 hotspots per Earth area; Lower 
mantle P-velocity (Clayton and Comer, 1983), 10 4 km/sec. Geopotentials may 
be converted to geoid elevations by dividing by the gravitational acceleration 
at the surface. 

Figure 3 - Degree-by-degree correlations, r t , between the hotspot distribution and 
the slab residual geoid (solid line). Correlations between the slab geoid model 
(Hager, 1984) and the observed geoid are shown for comparison (dashed line). 
Contours give the confidence of correlation, with a confidence level of 0.98 
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indicating only a 2% chance that the correlation is random. 

Figure 4 - Degree 2 comparison of the slab residual geoid, the vertically averaged 
lower mantle P-wave velocity model (Clayton and Comer, 1983), and the 
hotspot density distribution. Geoid lows, slow velocity anomalies, and low 
hotspot density areas are shaded. Hotspot density contour intervals are in 
units of 16.5 hotspots per Earth area. 

Figure 5 - Degree 6 comparison of the slab residual geoid, upper mantle shear 
velocity (200-500 km: Woodhouse and Dziewonski, 1984), and the hotspot den- 
sity distribution. Geoid lows, fast velocity anomalies, and low hotspot density 
areas are shaded. Hotspot density contour intervals are in units of 33 hotspots 
per Earth area. 

Figure 6 - Best-fitting hotspot/geoid response curve (solid dots) according to equa- 
tion 3 of text. Error bars (2a) indicate the uncertainty of the fit. Open circles 
give the best-fitting curve under the assumption of equal noise in both the 
slab residual geoid and the hotspot distribution (error bars not shown). 

Figure 7 - Slab residual geoid in harmonic degree bands 4r6, 7-12, and 10-20. 
Hotspot locations are shown by black dots. 

Figure 8 - Dynamic response functions, G t (r), for surface density contrasts of 
spherical harmonic degrees 2, 4, 6, and 8 plotted against radius, r , for six 
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Earth models. Models U, for uniform composition, left, permit mantle-wide 
flow; models C, right, have a (chemical) barrier at 670 km depth, causing 
stratification into separate upper and lower mantle flow systems. Models Ul 
and Cl have uniform viscosity; models UlO and CIO have a factor of 10 
viscosity increase below 670 km; models U100 and C100 have a factor of 100 
increase. The geoid responses are normalized to the geoid which would be 
obtained if the harmonic density contrasts were placed at the top surface with 
no dynamic compensation allowed. 

Figure 9 - Residual depth anomaly map of the Hawaiian swell (adapted from 
Schroeder, 1984). Flow stagnation contours for 300°K and 1000°K kinematic 
plume models are shown by bold lines. Cross sections AA’ and BB’ are plot- 
ted in Figure 20. 

Figure 10 - Example finite element grid used in plume and slab models. 

Figure 11 - Illustration of the A, B, and C plume models. Horizontal and vertical 
scales are equal. r) 0 is the reference viscosity. 

Figure 12 - Pseudo geoid response functions (see text) for the first three harmonics 
(wavelengths corresponding to degrees 4,8,12) of the cylinder (or box) for the 
finite element model A plume (solid line) and the high viscosity slab (dashed 
line). The dotted line gives the analytical solution for uniform mantle viscos- 
ity. 
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Figure 13 - Excess geoid signal perturbation relative to the isoviscous case, for the 
range of excess plume temperatures and viscosity exponents, 0 , tested at 
plume radius r 0 = 70 km. The 1 = 4 , 8, and 12 curves are for the first three 
harmonics of the cylinder for halo model A. 

Figure 14 - Buoyancy flux (r 0 =70 km) as a function of excess plume temperature 
and 0 for halo model A. The background viscosity is 10 21 Pa-sec 

Figure 15 - Buoyancy flux as a function of plume radius for halo model A. The 
right-hand vertical scale is linear and gives the fourth root of the flux. Solid 
lines are for excess temperatures of 400°K and 800°K with 0 — 28 . Dashed lines 
are for the isoviscous models ( 0 = 0 ). 

Figure 16 - Geoid signal perturbation (excess), relative to the isoviscous geoid, as a 
function of radius and viscosity contrast. The 1 = 4 , 8, and 12 curves are for 
the first three harmonics of the cylinder for the pipe model. 

Figure 17 - Buoyancy flux as a function of radius and viscosity contrast for the 
pipe model. The flux is normalized to the square of the excess temperature for 
density p 0 =5.14 and volume coefficient of thermal expansion 3z 10 S /°K. The 
right-hand scale is linear and gives the fourth root of the flux. Dashed lines 
give the theoretical flux for an infinite rigid pipe (mantle) for the same interior 
pipe viscosities and buoyancy forces (see text). 
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Figure 18 - Comparison of global hotspot response curve (solid line with error 
bars) and several dynamic models. Solid triangles give the / =4,8,12 geoid sig- 
nal for plume model B (/?=35). Solid diamonds are for model C (/?= 35). 
Model C ' (ft—O), indicated by open diamonds, has no viscosity contrast 
between the mantle and the plume. The “spherical Earth” model (dashed 
lines) was calculated analytically using the method of Richards and Hager 
(1984) with a 200 km radius plume, average temperature contrast 300°K, and 
no viscosity contrast. The dotted line gives the (spherical) analytical calcula- 
tion for the geoid spectrum from lithospheric delamination. 

Figure 19 - Pseudo geoid response functions (see text) for the model C plume 
(solid line), subducted slab (dashed line), and slab with weakened lithosphere 
(dash-dot line). The dotted line gives the analytical solution for the purely lay- 
ered response. 

Figure 20 - Comparison of dynamic topography from Models B and C with the 
observed Hawaiian swell topography. The heavy dotted line is from section 
A-A ' of Figure 9, and the lighter dotted line is from section B-B ' 500 km 
WNW of Hawaii (“anomalous” bathymetry adapted from Schroeder, 1984). 
Light and heavy solid curves are for strongly temperature dependent viscosity 
(/?= 35) in models B and C, respectively. Light and heavy dashed curves are 
for no viscosity contrast between plume and mantle (/3= 0) in models B and C, 
respectively. Deformation profiles for the theoretical models decay toward the 
zero level out to 5000 km from the plume axis, but are truncated in this figure 
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in order to emphasize deformation near the swell itself. 

Figure Bl - Illustration of the geometry for the high viscosity slab calculations. 
Horizontal and vertical scales are equal. 

Figure B2 - Pseudo response curves (/ =5,10) for models LlOO.u (dotted line), 
LlOO.slab (solid line), Ll00.slab200 (long dashes), and Ll00.slab400 (short 
dashes). 

Figure B3 - Pseudo response curves (1=5,10) for models LUlOO.u (dotted line), 
LUlOO.slab (solid line), LUl00.slab200 (long dashes), and LUlO0.slab40O (short 
dashes). 
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Figure 3 
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Residual Geoid: degree 2 



Lower Mantle P-Wave Velocity: degree 2 
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Figure 4 
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Residual Geoid: degree 6 



^ 0 j contour interval: 5 m 

Upper Mantle Shear Velocity: degree 6 
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Hotspot Density: degree 6 
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Figure 5 






- 76 - 



Figure 6 
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Residual Geoid: degree 4-6 



Residual Geoid: degree 7-12 



Residual Geoid: degree 10-20 



contour interval: 5 m (c) 


Figure 7 
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Figure 10 
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